%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file simulates the small model

clc; clear all; close all;

UT=1; % 0=production damage, 1=utility damage

%% Laissez fare allocation (Production damage)
Sbar=600;            % pre-industrial pollution level
S=900;               % pollution in excess of Sbar
nu=0.04;             % elasticity of output to energy
phi=0.48;            % elasticity of pollution to emission
gamma=5.7*10^(-5);   % production externality or linerar utility externality
gamma1=0;             % linear utility externality
%gamma2=0;
gamma2=2*gamma/S;   % quadratic utility externality to have same damage under LF as the production damage
ni=0.1;             % relative size of domestic economy
T=3*log((S+Sbar)/Sbar)/log(2);
E_lf=S/phi;          % energy under laissez faire
D_lf=1-exp(-gamma*S);% damage under LF
Y_lf=(1-D_lf)*E_lf^(nu); % output under LF
zeta=nu*Y_lf/E_lf;   % extraction cost
tauz=0;              % foreign tax

usedfont='Georgia';
FigW=35;
FigH=FigW;



%% Competitive equilibrium given optimal tax
x0=[E_lf,E_lf];
options = optimoptions('fsolve','MaxFunEvals',300000,'MaxIter',30000,'TolFun',1e-15,'Display','none','Algorithm','levenberg-marquardt');


[x,~,flag] = fsolve(@(xx) find_opt_tax(xx,gamma,gamma1,gamma2,phi,1,nu,zeta,tauz,UT,0),x0,options);
if flag<1
disp(['Computation fails']);
return
end
Ei=x(1);
Ez=x(2);


if UT==0
D=1-exp(-gamma*phi*(ni*Ei+(1-ni)*Ez));
Yi=(1-D)*Ei^(nu);
Yz=(1-D)*Ez^(nu);
taui=gamma*phi*ni*Yi;


else
Yi=(1-0)*Ei^(nu);
Yz=(1-0)*Ez^(nu);
taui=(Yi-zeta*Ei)*phi*ni*(gamma1+gamma2*(ni*Ei+(1-ni)*Ez));
end


%% Nash equilibrium
% Step 1: compute the optimal response of the domestic economy to the tax set by the rest of the world
taus=(0:1.0000e-8:1.0000e-04)';
N=length(taus);
tauis=zeros(N,1);
tauzs=zeros(N,1);


for j=1:N
[x,~,flag] = fsolve(@(xx) find_opt_tax(xx,gamma,gamma1,gamma2,phi,ni,nu,zeta,taus(j),UT,0),x0,options);
Ei=x(1);
Ez=x(2);


if UT==0
D=1-exp(-gamma*phi*(ni*Ei+(1-ni)*Ez));
Yi=(1-D)*Ei^(nu);
Yz=(1-D)*Ez^(nu);
tauis(j)=gamma*phi*ni*Yi;
else
Yi=(1-0)*Ei^(nu);
Yz=(1-0)*Ez^(nu);
tauis(j)=(Yi-zeta*Ei)*phi*ni*(gamma1+gamma2*(ni*Ei+(1-ni)*Ez));
end

if flag<1
disp(['Computation fails']);
return
end

% Step 2: compute the optimal response of the foreign economy to the tax set by the domestic economy
% It is like we assume that country i has size (1-ni), to compute the optimal tax of the other country
ni=1-ni;
[x,~,flag] = fsolve(@(xx) find_opt_tax(xx,gamma,gamma1,gamma2,phi,ni,nu,zeta,taus(j),UT,0),x0,options);

Ei=x(1);
Ez=x(2);


if UT==0
D=1-exp(-gamma*phi*((ni*Ei+(1-ni)*Ez)));
Yi=(1-D)*Ei^(nu);
Yz=(1-D)*Ez^(nu);
tauzs(j)=gamma*phi*ni*Yi;
else
Yi=(1-0)*Ei^(nu);
Yz=(1-0)*Ez^(nu);
tauzs(j)=(Yi-zeta*Ei)*phi*ni*(gamma1+gamma2*(ni*Ei+(1-ni)*Ez));
end

if flag<1
disp(['Computation fails']);
return
end
% Reset ni to the right value
ni=1-ni;

end

%Step 3 taus and tauzs are the reaction functions: approximate them to a grid of carbon taxes taus
tauis_app = interp1(taus,taus,tauis,'nearest');
tauzs_app= interp1(taus,taus,tauzs,'nearest');
reaction=[tauis_app, tauzs_app, taus];

%Step 4: find the Nash equilibrium starting from ne1 and ne2
ne1=reaction(1,1);
ne2=reaction(1,2);
A=1000;
j=1;
while A>0
ind1=find(taus==ne2(j));   % Take this action of foreign economy that constitutes a NE
ne1(j+1)=reaction(ind1,1); % this is the best response of Home...
ind2=find(taus==ne1(j+1)); 
ne2(j+1)=reaction(ind2,2); % and this is the best response of Foreign to the best response of Home
A1=abs(ne1(j+1)-ne1(j));   
A2=abs(ne2(j+1)-ne2(j));
A=norm([A1;A2]);           % Do we reach convergence? If yes, we found a NE, if not, go on with next action
j=j+1;
end
ne=[ne1(end),ne2(end)];

%% Coordination
ni=1;
[x,~,flag] = fsolve(@(xx) find_opt_tax(xx,gamma,gamma1,gamma2,phi,ni,nu,zeta,tauz,UT,0),x0,options);
Ei=x(1);
Ez=x(2);


if UT==0
D=1-exp(-gamma*phi*(ni*Ei+(1-ni)*Ez));
Yi=(1-D)*Ei^(nu);
Yz=(1-D)*Ez^(nu);
tau_cord_temp=gamma*phi*ni*Yi;
tau_cord=interp1(taus,taus,tau_cord_temp,'nearest');
tau_star=tau_cord;
save tau_star tau_star
else
Yi=(1-0)*Ei^(nu);
Yz=(1-0)*Ez^(nu);
tau_cord_temp=(Yi-zeta*Ei)*phi*ni*(gamma1+gamma2*(ni*Ei+(1-ni)*Ez));
tau_cord=interp1(taus,taus,tau_cord_temp,'nearest');
end
load tau_star
disp(['-----------------------------------------------------------------------']);
disp(['Nash equilibrium:         Home Tax ',num2str(round(ne(1)/tau_star,3)),', Foreign Tax ',num2str(round(ne(2)/tau_star,3))]);
disp(['Coordination equilibrium: Home Tax ',num2str(round(tau_cord/tau_star,3)),', Foreign Tax ',num2str(round(tau_cord/tau_star,3))]);
disp(['-----------------------------------------------------------------------'])

% disp(['-----------------------------------------------------------------------']);
% disp(['Nash equilibrium:         Home Tax ',num2str(round(ne(1)*100,7)),'%, Foreign Tax ',num2str(round(ne(2)*100,7)),'%']);
% disp(['Coordination equilibrium: Home Tax ',num2str(round(tau_cord*100,7)),'%, Foreign Tax ',num2str(round(tau_cord*100,7)),'%']);
% disp(['-----------------------------------------------------------------------'])



%% Plotting
if UT==0
tauis_prod=tauis;
tauzs_prod=tauzs;
save tau_prod tauis_prod tauzs_prod
end

if UT==1 && gamma1==gamma
tauis_lin=tauis;
tauzs_lin=tauzs;
save tau_lin tauis_lin tauzs_lin
end

if UT==1 && gamma2==gamma/S
tauis_quad=tauis;
tauzs_quad=tauzs;
save tau_quad tauis_quad tauzs_quad
end

if UT==1 && gamma2==2*gamma/S
load tau_prod
load tau_lin
load tau_quad   
usedfont='Georgia';
FigW=20;
FigH=20;

figure;
set(gcf,'color', 'white',...
        'PaperUnits','centimeters','PaperSize',[FigW FigH],...
        'PaperPosition',[0,0,FigW,FigH],'PaperPositionMode','manual',...
        'Units','centimeters',...
        'Position',[0,0,FigW,FigH]);
subplot(2,2,1)
plot(taus/tau_star,tauis_prod/tau_star,'Linewidth',2)
axis tight
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
title('Production externality','Fontsize',14,'fontweight','bold')
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')

subplot(2,2,2)
plot(taus/tau_star,tauis_lin/tau_star,'Linewidth',2)
%axis tight
axis([taus(1)/tau_star taus(end)/tau_star tauis_lin(1)/tau_star-0.001, tauis_lin(1)/tau_star+0.001])    
title('Linear disutility','Fontsize',14,'fontweight','bold')
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
subplot(2,2,3)
plot(taus/tau_star,tauis_quad/tau_star,'Linewidth',2)
axis tight
title('Quadratic disutility, low \gamma_2','Fontsize',14,'fontweight','bold')
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
subplot(2,2,4)
plot(taus/tau_star,tauis/tau_star,'Linewidth',2)
axis tight
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
title('Quadratic disutility, high \gamma_2','Fontsize',14,'fontweight','bold')
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
print('-dpng','-r500',['./figures/reaction.png'])

return

figure;
set(gcf,'color', 'white',...
        'PaperUnits','centimeters','PaperSize',[30, 10],...
        'PaperPosition',[0,0,30, 10],'PaperPositionMode','manual',...
        'Units','centimeters',...
        'Position',[0,0,30, 10]);
plot(taus/tau_star,tauis_prod/tau_star,'Linewidth',2)
hold on
plot(taus/tau_star,tauis_lin/tau_star,'r:','Linewidth',2)
hold on
plot(taus/tau_star,tauis_quad/tau_star,'g--','Linewidth',2)
axis tight
plot(taus/tau_star,tauis/tau_star,'k--','Linewidth',2)
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
axis tight
legend('Production externality','Linear disutility','Quadratic disutility, low \gamma_2','Quadratic disutility, high \gamma_2','Fontsize',14,'location','northeast')
print('-dpng','-r500',['./figures/reaction_compact.png'])


figure;
set(gcf,'color', 'white',...
        'PaperUnits','centimeters','PaperSize',[FigW FigH],...
        'PaperPosition',[0,0,FigW,FigH],'PaperPositionMode','manual',...
        'Units','centimeters',...
        'Position',[0,0,FigW,FigH]);
subplot(2,2,1)
plot(taus/tau_star,tauzs_prod/tau_star,'Linewidth',2)
axis tight
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
title('Production externality','Fontsize',14,'fontweight','bold')
subplot(2,2,2)
plot(taus/tau_star,tauzs_lin/tau_star,'Linewidth',2)
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
%axis tight
axis([taus(1)/tau_star taus(end)/tau_star tauzs_lin(1)/tau_star-0.001, tauzs_lin(1)/tau_star+0.001])    
title('Linear disutility','Fontsize',14,'fontweight','bold')
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
subplot(2,2,3)
plot(taus/tau_star,tauzs_quad/tau_star,'Linewidth',2)
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
axis tight
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
title('Quadratic disutility, low \gamma_2','Fontsize',14,'fontweight','bold')
subplot(2,2,4)
plot(taus/tau_star,tauzs/tau_star,'Linewidth',2)
ax = gca;
ax.FontSize = 14; 
ax.YAxis.Exponent = 0;
axis tight
ax.YAxis.Exponent = 0;
xlabel('\tau_F','Fontsize',14,'fontweight','bold')
ylabel('Optimal \tau_H','Fontsize',14,'fontweight','bold')
title('Quadratic disutility, high \gamma_2','Fontsize',14,'fontweight','bold')






end